##########################################
###Autocratic Elections
###Stabilizing Tool or Force for Change
###Replication of GAM models
#########################################
rm(list = ls())
library(foreign);library(pspline);library(survival);library(MCMCpack)
library(ZeligGAM);library(R2jags)
library(mgcv);library(glmnet)
library(bootstrap);library(plyr);library(stargazer)

data <- read.dta("ReplicationData.dta")
dataallelections <- read.dta("ReplicationDataAllElections.dta")


tempdata <- subset(data,select=c(gwf_fail,timesinceelection, timesincelegislativeelection, timesinceexecutiveelection, timesincecontestedelection,timesinceuncontestedelection, regiondem, loggdp, gdp_grow , 
                                   cow_milsize , resdep2,  duration, duration2, duration3,region,decade))
tempdata <- na.omit(tempdata)
gam.model <- gam(gwf_fail ~ s(timesinceelection, k=4,fx=TRUE,bs="tp")  + regiondem + loggdp + gdp_grow + 
                   cow_milsize + resdep2 +  duration + duration2 + duration3,
                     family=binomial(link = "logit"), data=tempdata)
summary(gam.model)
stargazer(gam.model, 
          type="latex", single.row=TRUE,
          out="../Results/GAM1Results1.tex")

pdf("../Figures/GAMeffects1.pdf")
plot(gam.model, rug=TRUE, se=TRUE,
     select=1, main="",ylab="Partial effect of election",
     xlab="Time since election", shade=TRUE, xlim=c(0,20))
abline(0,0)
dev.off()

for(i in 4:8){
  gam.model <- gam(gwf_fail ~ s(timesinceelection, k=i,fx=TRUE,bs="tp")  + regiondem + loggdp + gdp_grow + 
                     cow_milsize + resdep2 +  duration + duration2 + duration3,
                   family=binomial(link = "logit"), data=tempdata)
  name <- paste("../Figures/GAMeffectsDF",i,".pdf",sep="")
  title <- paste("Degrees of freedom=",i,sep="")
  pdf(name)
  plot(gam.model, rug=TRUE, se=TRUE,
       select=1, main=title,ylab="Partial effect of election",
       xlab="Time since election", shade=TRUE, xlim=c(0,20))
  abline(0,0)
  dev.off()
}


##Figure 1 on all types of elections

tempdataall <- subset(dataallelections,select=c(gwf_fail,timesinceelection, timesincelegislativeelection, timesinceexecutiveelection, timesincecontestedelection,timesinceuncontestedelection, regiondem, loggdp, gdp_grow , 
                                 cow_milsize , resdep2,  duration, duration2, duration3,region,decade))
tempdataall <- na.omit(tempdataall)
gam.model <- gam(gwf_fail ~ s(timesinceelection, k=4,fx=TRUE,bs="tp")  + regiondem + loggdp + gdp_grow + 
                   cow_milsize + resdep2 +  duration + duration2 + duration3,
                 family=binomial(link = "logit"), data=tempdataall)

summary(gam.model)

stargazer(gam.model, 
          type="latex", single.row=TRUE,
          out="../Results/GAM1ResultsAllElections.tex")

pdf("../Figures/GAMeffectsAllElections.pdf")
plot(gam.model, rug=TRUE, se=TRUE,
     select=1, main="",ylab="Partial effect of election",
     xlab="Time since election", shade=TRUE, xlim=c(0,20))
abline(0,0)
dev.off()

##with decade and region dummies
gam.model.2 <- gam(gwf_fail ~ s(timesinceelection, k=4,fx=TRUE,bs="tp")  +  regiondem + loggdp + gdp_grow + 
                   cow_milsize + resdep2 +  duration + duration2 + duration3 + as.factor(region) + as.factor(decade),
                 family=binomial(link = "logit"), data=tempdata)
summary(gam.model.2)
pdf("../Figures/GAMeffectsRegionDecade.pdf")
plot(gam.model.2, rug=TRUE, se=TRUE,
     select=1, main="",ylab="Partial effect of election",
     xlab="Time since election", shade=TRUE, xlim=c(0,20))
abline(0,0)
dev.off()

##UNCONTESTED AND CONTESTED ELECTIONS 
gam.model.contested <- gam(gwf_fail ~ s(timesincecontestedelection, k=4,fx=TRUE,bs="tp")  + regiondem + loggdp + gdp_grow + 
                   cow_milsize + resdep2 +  duration + duration2 + duration3,
                 family=binomial(link = "logit"), data=tempdata)
summary(gam.model.contested)
gam.model.uncontested <- gam(gwf_fail ~ s(timesinceuncontestedelection, k=4,fx=TRUE,bs="tp")  + regiondem + loggdp + gdp_grow + 
                   cow_milsize + resdep2 +  duration + duration2 + duration3,
                 family=binomial(link = "logit"), data=tempdata)
summary(gam.model.uncontested)
pdf("../Figures/GAMeffectsContestedvsUncontested.pdf")
par(mfrow=c(2,1))
plot(gam.model.contested, rug=TRUE, se=TRUE,
     select=1, main="Contested elections",ylab="Partial effect of election",
     xlab="Time since election", shade=TRUE, xlim=c(0,20))
abline(0,0)
plot(gam.model.uncontested, rug=TRUE, se=TRUE,
     select=1, main="Uncontested elections",ylab="Partial effect of election",
     xlab="Time since election", shade=TRUE, xlim=c(0,20))
abline(0,0)
dev.off()

##EXECUTIVE AND LEGISLATIVE ELECTIONS
gam.model.executive <- gam(gwf_fail ~ s(timesinceexecutiveelection, k=4,fx=TRUE,bs="tp")  + regiondem + loggdp + gdp_grow + 
                             cow_milsize + resdep2 +  duration + duration2 + duration3,
                           family=binomial(link = "logit"), data=tempdata)
summary(gam.model.executive)
gam.model.legislative <- gam(gwf_fail ~ s(timesincelegislativeelection, k=4,fx=TRUE,bs="tp")  + regiondem + loggdp + gdp_grow + 
                               cow_milsize + resdep2 +  duration + duration2 + duration3,
                             family=binomial(link = "logit"), data=tempdata)
summary(gam.model.legislative)
pdf("../Figures/GAMeffectsExecutiveLegislative.pdf")
par(mfrow=c(2,1))
plot(gam.model.executive, rug=TRUE, se=TRUE,
     select=1, main="Executive elections",ylab="Partial effect of election",
     xlab="Time since election", shade=TRUE, xlim=c(0,20))
abline(0,0)
plot(gam.model.legislative, rug=TRUE, se=TRUE,
     select=1, main="Legislative elections",ylab="Partial effect of election",
     xlab="Time since election", shade=TRUE, xlim=c(0,20))
abline(0,0)
dev.off()


